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ANALYSIS OF THE TRANSIENT 


BEHAVIOR OF RUBBING COMPONENTS 


ABSTRACT 


Finite element equations are developed for studying 
deformations and temperatures resulting from frictional 
heating in sliding system. The formulation is done for 
linear steady state motion in two dimensions. The 
equations include the effect of the velocity on the 
moving components. This gives spurious oscillations in 
their solutions by Galerkin finite element methods. A 
method called "streamline upwind scheme" is used to try 
to deal with this deficiency. The finite element 
program is then used to investigate the friction of 
heating in gas path seal. 
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LIST OF SYMBOLS 


0^ Specific heat at constant pressure 
E Young's modulus 

F Body force 

H Hilbert spaces 

k Thermal conductivity 

Artificial conductivity 
N Shape function 

q Heat flux 

T Absolute temperature 

U Displacement 

V Velocity 

W Weighting function 

a Thermal expansion coefficient 

A Lame's constant 

p Shear modulus 

p Density 

£,r) Natural coordinates 

v Poisson's ratio 

uj Angular velocity 
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Chapter 1 
INTRODUCTION 

The first law of thermodynamics expresses the energy balance 
during mechanical and thermal process. In the analysis of rubbing 
problem, the loss of mechanical energy (frictional energy) is 
transformed in its largest percentage to thermal energy. During 
high speed sliding, contact patches are formed. An analytical 
treatment of stresses (or displacement) and temperature distribution 
near the contact patches is necessary. A transient finite element 
heat conduction analysis has shown (ref. 14) that within a very short 
time after establishment of the contact zone the temperature 
distribution approached a steady state relative to a stationary 
observer. The length of time required to reach this quasi-steady 
state is so short that it may be concluded that within the contact 
patches a steady state temperature distribution occurs. Therefore 
it is not necessary to do a transient temperature analysis. A finite 
element formulation will be done for linear steady state. The 
formulation will be given in a form that could be expanded to 
inelastic, non-linear problems. 

The advantage of the Finite Element Method is that it is 
possible to model finite geometry of complex shapes or different 
material properties. Both temperature and stress analysis could be 
done by similar modeling. Indeed, a thermomechanical analysis could 
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be carried out using one element grid and two linked finite element 
programs. The major difficulty in applying a Finite Element Method 

is that the convection operators are nonsymmetric . For instance the 
Galerkin Finite Element Method is successful when applied to linear 
symmetric operators, but these methods usually give spurious 
oscillations in their solutions when applied to convection dominated 
problems. A "streamline upwind scheme" (ref. 10) is used to deal 
with this problem by adding an artificial conductivity in a manner 
which stabilize the solution without destroying the physics of the 
problem. 

In this work, a review of literature about the principal 
subjects is given, followed by a formulation of the weak form for the 
heat transfer equation and the thermoelasticity equation. Then a 
finite element formulation is developped for both thermal and 
thermoelastic equation for a two dimension solid. The resulting 
finite element program, which gives the displacement and the 
temperature distribution, is first compared to an analytical solution 
such as a semi-infinite plane under a heat flux (ref. 1). The 
program is then used to a problem of rubbing contact at high 
velocities in a gas path seal. 


Chapter 2 
LITERATURE REVIEW 

The heat transfer theory started with Fourier's law of heat 
conduction: 

, a dT 

1 - ' k A ST 

When the body is moving with a given velocity, a convection 
term is added to the equation (ref. 13 & 17). The use of Galerkin 
Finite Element Method to solve the heat problems for a moving body 
give rise to spurious oscillations. These oscillations can be 
removed in this case by severe mesh refinement which undermines the 
practical utility of the methods (ref. 8, 9 and 10). 

New schemes were developped trying to deal with this deficiency. 
The first scheme appeared by Roache in ref. 6 as a classical upwind 
difference scheme. It has been noted that the Galerkin Finite 
Element Method produces central difference type approximations to the 
advection (conduction) term. In finite difference theory, the 
adverse behavior of central differences in these circumstances has 
long been noted. But this method was considered as inaccurate. 

Heinrich proposed a new scheme (ref. 8). The Finite Element 
Method is applied using weighted residual formulation with 
bilinear quadrilateral element shape functions, and non-symmetric 
weighting functions which are different from the shape 
functions, and depend on parameters which allow the amount of 
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"upwinding" to be controlled. An increase in accuracy could be 
obtained by varying these parameters from element to element. This 
method is now known as Petro-Galerkin method. But the two 
dimensionnal quadrilaterals proposed by ref. 8 distord the diffusion 
(conduction) term when upwinding is applied. It seems very difficult 
to find an upwinding function that does not disturb the diffusion 
(conduction) operator, yet upwinds the advection (convection) term. 
Hansen and Von Flotow (ref. 16) noted that it might be better to 
apply upwind weighting to the advection (conduction) term only and 
central to the remainder of the equation. 

Another simpler technique proposed by Hushes (ref. 9) called 
quadrature upwinding was based on moving the integration points in 
the Galerkin Finite Element Method. But he came later with Brook 
(ref. 10) to propose a new multi-dimensionnal upwind scheme. The 
method was applied successfully to one dimension, and then 
generalized for two dimensions. This method, called "streamline 
upwind scheme", is applied to the advection-dif fusion equation, and 
then to Navier-Stokes equations. 

The thermoelasticity equation is derived from the principles of 
thermodynamics. A simple formulation of the equation is presented 
in ref. 2 & 3 as: 

yU. . . + (A+y) U. - (3A+2y) aT . + F, = 0 
i.33 3.31 .i i 

Together with the heat transfer equation, the thermoelasticity 
equation (two equations in dimension two) leads to the determination 
of the displacement and temperature fields in a thermomechanical 
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problem such as the problem of rubbing contact at high sliding 
velocities . 

The thermal analysis of bodies in sliding contact has attracted 
the interest of many investigators because of its importance in many 
situations in which friction occurs: bearing, seals, brakes, 
clutches,... Various methods have been proposed, but none has proven 
universally acceptable. Many surface temperature analyses have been 
based on heat source methods (ref. 1), in which the solution for 
temperature distribution due to a point source on a surface is used 
to develop the solution for a distributed heat flux within a contact 
patch on the surface of an infinite half space. The difficulties 
involved with application of heat source techniques to bodies of 
finite dimensions led to the development of integral transform 
technique presented by Ling (ref. 4). Although this method have 
been successfully applied to a number of problems with different 
geometries, its limitations to simple shapes and its mathematical 
sophistication have kept this from being widely used by engineers. 

Kennedy tried to solve the problem of rubbing contact at high 
sliding velocities by using the Finite Element Method. He applied 
the problem to two examples: aircraft disk brakes and gas path seals 
in turbine engine (ref. 11). The first examples was presented 
before in the study of transient temperature in disk brakes in 
ref. 5, considered as one of the first documention that use Finite 
Element Method in such problems. He retreated the second example 
experimentally and analytically in ref. 14 & 15. In ref. 11, 
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Kennedy used one element grid and two linked finite element programs 
to make a thermomechanical analysis of the contact. 


Chapter 3 

FORMULATION OF THE PARTIAL 
DIFFERENTIAL EQUATION 


3-1. Heat transfer: 

Let k , the thermal conductivity, be constant, let P, the 
density, and C^, the specific heat be constant. Let q be the heat 
flux. In a steady state, the heat transfer partial differential 
equation for a moving body with a constant velocity V. is: 

k.. T ..-pC V. T . + q = 0 
ij ,ij P i .i H 

The first term represents the conduction heat transfer. The 
second term represents the convection heat transfer. The problem 
defined in the equation above could be given after applying 
Galerkin's method, as: 

Find T G H 2 for all W G H° scuh that: 


[ [ W k. . T . . - W p C V. T.+Wq]dft = 0 

J Q .ij P i .i 

Where T is assumed to satisfy the essential boundary 
conditions. An integration by parts allows us to rewrite the problem 
as: 

Find T G H 1 for all W G H 1 such that: 



W 


k. .. 
ij 


+ W P C 


V. 

l 


T . •] dfl 


f W q dfl 
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Here homogenous natural boundary conditions have been assumed 
in those sections of the boundary where T is not specified. 

3-2. Thermoelasticity: 

The constants A and P are respectively Lame's constant and shear 
modulus. They are related to Young's modulus E and Poisson's ratio 
v by: 

v E E 

A = *, P = 

(1 + v) (1 - 2 v) 2 (1 + V) 

The Navier's equation, with temperature changes, in terms of 

displacement U is given by: 

p U. . . + (X + p) U . .. - (3A + 2 p) aT . + F, = 0 
i.JJ ' J.Ji i 

Where a is the coeficient of the thermal expansion and is the body 

force. 

Applying the Galerkin method to the equation above within an 
element results in the following formulation: 

Find U 6 H 2 , T 6 H 1 for all W6H° such that: 

f [ W p U. . . + W (X + p) U. .. - W (3 A + 2 p ) a T . 
Jq *OJ JO 1 » x 

+ W F. ] dS7 = 0 

l J 

After an integration by parts, we write the formulation as: 

Find U6H 1 , T6H 1 for all W6 h’ such that: 

( [W . p l). ,+W . ( X + P ) U . . + W (3X + 2p) a T,i ] dfi 
1 ,J ifJ .i JO 

= [ W F. dfi 
'ft 1 
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This equation, called the weak form, could be obtained if the 
law of conservation of energy is used. 


Chapter 4 

FINITE ELEMENT FORMULATION 

4-1. Heat transfer: 

4-1-1. Streamline upwind scheme: 

The heat transfer equation has a convection operator which has 
been shown to carry spurious oscillations in finite element solution 
(ref. 10). The conventional approach to mitigate these oscillations 
is to introduce an artificial diffusion (conduction) term in the 
heat transfer equation. This method is called "streamline upwind 
scheme". 

The weak form becomes: 


[ [ W . (k . . 4k..) T .+WpC V.T. ]dn=f W q dfi 

Jfi fl 1J 3 P i .1 

A 

If the artificial conductivity k„ is correctly chosen, no 
oscillation will occur in the Galerkin Finite Element formulation. 
In ref. 10, the following technique is presented: 


Assume k . . = k u . u . 

iJ i J 


where u . = with Mult 2 = u . u . and u . = p C V . 

i IL II ii i P i 


A 

k is a scalar artificial conductivity. 

Assume that the coordinates are chosen such that locally x x - 
direction is aligned with streamlines and the x 2 -direction is 
perpendicular. Then the artificial conductivity matrix in this 
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coordinate system is: 



In case of bilinear quadrilateral, in two dimension, k is chosen 

as: 

k = i ( £' u c + n ) 

A A 

£ andndefine the location of the quadratic point and are given by: 




Figure 1. Typlcel four-sode quedriUterel finite ele- 
ment geoaetry. 
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u and k are evaluated at the origin of the element in£~n 
coordinates, could be expressed by: 


a . = 
1 


PC V. h. 
p 1 1 


; i = c , n 


4-1-2. Heat transfer finite element formulation: 

The weak form is know given by: 

[ [W. (k. . + k. .) T , + WnC V.T. ] dl]= [ W q dfi 

If the body is divided into a number of finite elements, the 
weighting function W and the temperature distribution T within a 
discrete element may be approximated by: 

W - Wj Nj 

T * T i h 

Where Nj is the shape (interpolation) function of the element. 
Wj and Tj are constant at the nodes. The majuscule subscripts I 
indicate the node number. 

Substituting these approximations into the weak form we get: 

J l Wj htl Ck. . + i^) Tj N Jf . + Wj Nj (p Cp V.) Tj N Jtl ) dt> 

= f Wj Nj q dfl 
Jfi 

which may be simplified as follows: 

1 { ! N i,i < k ij + k ij> N j,j + N i «> C P V N J,i 1 dn,T J 

= [ N, q d ft 
Jo 1 
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[ Q ] = f Lqdfi 

These equations are the finite element equations which will be 
formed, assembled, and solved for the temperature distribution 
through the body. Because of non-symmetry of the second term 
in the stiffeness matrix [k], its presence require the use of 
solution routines different from those used in most finite element 
programs. 

To compare the solution obtained with the solution without the 

A 

use of the upwind scheme, we can just set = 0. The same finite 
element equations can also be used for both stationary and moving 
components of a sliding system by setting V = 0 for elements in the 
stationary body. 

4-2. Thermoelasticity: 

The weak form of the thermoelasticity equation was found to be: 

ft W ,j VI u 1 + w (X + V) U + W (3* + 2U ) a T ] dfl 
■'ft J - L, J j»J '- 1 
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W F. dft 
1 


The weighting function W, the displacement distribution U, and 
the temperature distribution T, within a discrete element may be 
approximated by: 

W = Wj Nj 


1) = Uj Nj 


T - Tj N, 

where the majuscule subscripts indicate the node number and N are 
the shape function defining the type of element. For some elements, 
the shape functions are given in the appendix. , Uj and are 

constants at the nodes. 

After substituting these approximations into the weak form, we 

get: 


I c f W X N I.J UU J N J,j + M I "l.i < A + U)U J N J.j 

+ W l“l < 3X + 2U) a Tj N Jf . } dn = } n Wj Nj Pj 
which may be simplified as follows: 




< J n lN i.j UN j,j + N i.i (A+,J)N j.j ldn,u j 

+ ! | I "l (3X + 2«)«* j?1 J d!5)Tj. | Nj F, 


dn 


'A - 1 - u 'n 

These equations can be written in the matrix notation as: 
[ K e 3 [ U ] + [ K et ] [ T ] = [ F ] 
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where [ ] = I [ N T . U N. . + N, . ( ^+U) N T . ] d£2 is the 

t- Jq ± , J J » J t , l J » J 


'n 

elastic stiffeness matrix. 

[ 0 ] = Uj 


[ k £T ] - 




[ N, (3^ + 21i ) Q N, . ] dfi is the coupled 

X *J t X 


stiffeness matrix. 
[ T ] = Tj 


[ F ] = 




Nj Fj dfi 


Based on the equations above and the heat transfer finite 
element equations, defined in the previous section, a finite element 
program has been written. This program, which will be used in a 
thermomechanical analysis, will give the displacement and the 
temperature distribution through the body. It should be noted that 
one element grid can be used. With this program simpler problems 
could be treated. For example: 

-Heat transfer without thermal stresses problem (an example of 
this problem is given later): Set no forces and a = 0. 

-Elasticity problem without temperature change: Set no flux. 

Some basic routines and a typical input data are given in the 


appendix . 


Chapter 5 

PROBLEMS AND SOLUTIONS 

5-1. Semi-infinite solid under heat friction: 

5-1-1. Description of the problem: 

The semi-infinite source strip is defined by: x£[-b,b); 
yG[0, °°]; in the plane z = 0. The heat is applied at the 
rate Q per unit time per unit area over the strip. The 
surrounding media moves across it with velocity V in the 
direction of the x-axis. This problem was treated by 
Carslaw and Jaeger (Ref. 1) using the heat source method. 


K-2b-^J __V 



Fig. 2-Semi-infinite solid under friction. 
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5-1-2. Finite element method: 

The problem defined above is treated by the finite 
element program. The elements used in the mesh are quadri- 
lateral (linear or quadratic). The heat flux is applied 
uniformly distributed on a width of 2b at the nodes. Four 
meshes were used, which are: 

- mesh 1: 200 elements, 231 nodes, four-node quadrilateral 

element, with upwinding (Fig. 3). 

- mesh 2: 200 elements, 231 nodes, four-node quadrilateral 

element, without upwinding. 

- mesh 3: 300 elements, 981 nodes, eight-node quadrilateral 

element, with upwinding. 

- mesh 4: 460 elements, 1493 nodes, eight-node quadrilateral 

element, with upwinding. 

5-1-3. Results: 

Carslaw and Jaeger (Ref. 1) gave the temperature as the 
following integral: 

T * Tw PI K ° + u! > 1/2 

J x-B 


where Ko(x) is the modified Bessel function of the second 
kind of order zero and X, Z, B are dimensionless quantities 
introduced as: 


y Vx 7 Vz R Vb . ... „ _k_ 
X = 2 k * z = 2K ’ B = 2K ’ th K p C 


Some values of the surface temperature of the solid 


are shown in Fig. 4 for B = 10. The curve represents 


r^o 
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£ifl.3-Linear quadrilateral mesh for a "semi-infinite solid under 


heat friction". 



R«f. I 

.quadrat 460 ok 
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(vTkV/2KQ) vs. (x/b) given by Ref . 1. The results by finite 
element method, for each mesh, is plotted on the same 
representation . 

5-1-4. Discussion: 

When the quadratic elements (8 nodes) are used the curve 
(7iTkV/2KQ) vs. (x/b) approaches the one given by the heat 
source method (Ref. 1). When a smaller number of linear 
elements (4 nodes) is used the curve is less accurate but 
has the same pattern. The oscillations in the solution are 
minimized when the upwinding is used. But it is still 
necessary to make a mesh refinement in order to get accurate 
results. 

5-2. Gas path seals: 

5-2-1. Description of the problem: 

The gas path seals are used in the turbine engines of 
modern aircraft to prevent the axial flow of working fluid 
(air) around rotating engine components. Reduction of the 
clearances between rotating and stationary component of these 
seals can decrease the consumption of the specific fuel and 
increase the effeciency of the engine. However, such 
reduction in clearance may result in occasional rubbing 
between the rotating and the stationary seal components as 
engine deflection occurs. These rubs, which occurs at very 
high sliding speeds, can cause high surface temperatures, 
excessive wear of the seal components, and possible damage 
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to the engine. The development of gas path seal designs 
have been retarted by an incomplete understanding of the 
temperatures, stresses, and deformations which occur during 
high speed seal rubs (Ref. 11). 

An attempt to get the temperature and deformation in this 
sliding system is done using the finite element analysis. 

A model has been developed to simulate the rubbing contact 
which occurs in gas path seals between the rotating knife 
edge and a stationary seal segment. The outer gas path seal, 
assumed to have a circular section with 5 layers (Fig. 5 & 6), 
is rotating at 20,000 r.p.m.. The friction force is taken as 
100 lb.. The material properties are given in Table 1. On 
the external surface the temperature is taken to be 70°F and 
the displacement is zero. On the internal surface the flux, 
considered concentrated at the lower point is taken to be: 
q = N P V <= 100 x (0.3) x (20,000 x 2 ttx 4.85) 

= 1.8284 x 10 7 BTU/min in 2 
where u is the friction coefficient. 

5-2-2. Finite element method: 

The thermomechanical program developed before is used 
to the rubbing contract problem. A finite element mesh 
based on quadratic quadrilateral (8 nodes) elements, shown in 
Fig. 7, is used in the development of the model. 

5-2-3. Results: 


3-1. Temperature: 



Table 1 


Material properties 


Layer 

1 

2 

3 

4 

5 

Int. radius (in.) 

4.85 

4.91 

4.94 

4.97 

5.00 

Ext. radius (in.) 

4.91 

4.94 

4.97 

5.00 

5.11 

Material 

100 YSZ 

85% YSZ 
15% CoCr 
A1 Y 

70% YSZ 
30% CoCr 
A1 Y 

40% YSZ 
60% CoCr 
A1 Y 

MAR-M- 

50g 

Young's modulus 
(lb/in 2 )xl0 6 

2.00 

2.00 

8.00 

17.75 

15.60 

Poisson's ratio 

0.25 

. 0.26 

0.27 

0.28 

0.30 

Coef. of expansion 
(in/ in°F)xlO 6 

4.83 

7.70 

8.38 

9.52 

12.20 

Thermal cond. . 

(BTU/min in°F)xlO 

8.11 

16.30 

20.95 

25.52 

422.98 

Density 

(lb/in 3 ) 

0.155 

0.180 

0.205 

0.254 

0.320 

Specific heat 5 
(BTU/lb°F)xl(j 

0.161 

0.161 

0.161 

0.158 

0.155 
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T=70°F, U=0. 



Fig 7- Finite element mesh for gas path seal friction problem. 
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The finite element results give the temperature at any 
node. A plot of the temperature on the internal surface vs. 
the angle near the contact is given in Fig. 8. A peak of 
temperature occurs at the location of the friction contact. 

3-2. Deformation: 

The finite element results give the displacement in 
x-direction and in y-direction of any node. The curve of the 
deformation in x-direction vs. the angle for the interior 
surface is given in Fig. 9. The displacement is positive for 
the nodes at the right of the friction contact and negative 
for the nodes at its left. The curve of the deformation in 
y-direction vs. the angle is given in Fig. 10. This curve 
presents a peak at the friction contact. 

5-2-4. Discussion: 

At the location of the friction contact, the temperature 
and the deformation are increased rapidly. The temperature 
might reach the melting point of the material. The 
temperature peak was predicted by Marsher (Ref. 12) and 
confirmed by Kennedy (Ref. 11). The latter conformed also 
that the maximum amount of deformation occurs on the contact 
surface, with magnitudes decreasing very rapidly in a 
direction normal to the surface. 

The oscillations that appear clearly in the x-direction 
displacement (midside node) are due to the presence of the 
oscillations in the temperature distribution. When the 
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Fig.8-Internal temperature distribution under friction near 
contact. 
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upwjnding term is increased, the magnitude of temperature 
and displacement decrease but the oscillations don't vanish. 


I 
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Chapter 6 
CONCLUSIONS 

A finite element program has been developed which can 
efficiently and accurately predict temperature and deformation 
distribution in a thermomechanical problem such as sliding systems. 
The method includes velocity effects with the use of the streamline 
upwind scheme which eliminates spurious oscillations and is therefore 
very useful in cases involving high speed sliding. The comparison 
with some results that used heat source method shows the accuracy 
of the method. A special application to gas path seal problem shows 
the existence of a peak in the temperature and deformation near the 
contact, indicating a possibility of melting and excessive wear 
which might provoke an early thermomechanical failure. Further 
studies could investigate the importance of thermal conductivity of 
the moving and stationary component on decreasing surface 
temperature and deformation. 
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Appendix I 
SHAPE FUNCTIONS Nj 


1-Linear quadrilateral element: 

N j = 4 (i - £ ) (i - n ) 

n 2 = i (l + £ ) (l -n ) 

N 3 - i (l +€ ) (l + n ) 

n 4 « i (i -5 ) (l +n ) 


♦ u 


4 (- 1 , 1 ) 


1 (-1,-1 \ 


— ® 3 ( 1 , 1 ) 


i2 (1,-1) 


2-Quadratic quadrilateral element: 

Nj = -t ( l- C )(i-n )(i+5 + n ) 
n 2 = -i (i+5 )( l- n )(i-5 + n ) 
n 3 = -i (i+5 )(i+n )d-C-n ) 

N 4 = -{ (1-5 )(l+n )(1+C-T1 ) 

n 5 = h (l-C 2 )(i- n ) 
n 6 = h d+5 )d-n 2 ) 
n 7 = i (i -£ 2 )(i+n ) 



n 8 = i (l-n 2 )(i- 5 ) 


ORIGINAL PAGE IS 
OF POOR QUALITY 


Appendix II 


STIFFENESS MATRIX ROUTINE 


OOOl 

COO.! 

0003 

OOOA 

0005 

CCOe 

0007 

0006 

000 ? 

00 jo 
con 
COS 2 
C 01 3 

001 A 
00 IS 
OOif 
0017 
OC 1 £ 
001 c 
00*0 
0031 
0033 
0033 

mt 

00 3 c 

002 ? 

m 

0030 

0031 
0033 
0033 
C 03 A 
0033 
003 * 
OC 3 7 
002t 
OC?A 
00 AO 
OCA 1 
00 A 3 
00 A 3 

mi 

OCA 6 
00 A 7 

wai 

0050 

POM 

com 


EUECCl/T 1«E Cf 1 IF C Jt . t , u r .Mor ir.l YM «£M 

imun AtAnt <a-k«o-z> 


005A 


lUM 

c*»»a 
cmi 
r#»» * 

C>>*4 

CI'M 

rim 

Ctm 

c»»*» 

Cun 

Cmi 

Cim 

dm 

cim 

t>m 

cmi 

C***i 

Cim 

tmi 

Cim 

rim 

(im 

Cim 

Cim 

cim 

Cim 

cmi 

cmm 

cim 

cmm 

cmi 


S u 6 k 0 u 1 


HI t fc T 1 f 


V' 

I 


this eurrouiinl c a l c ul ai r e tut valiif ur inr 

FOR ANT DM£ R TWO MPINFICIWAL 11 0FAKA«M M C 


IT Iff KC‘C 
Ft. l fit N i 


I'i A01 


* » I t 

im 

r.Aif.ix »»»» 

*» »t 

■ in 
im 
»» » • 
M »1 
«» *1 
im 
»»»t 
»» » » 
mi 
im 
mm 

i » • « 

mm 

• m 
im 
mi 

tut 
till 

mi 
mi 

t»» • 

mi 

• ait 
Mt* 
»*»» 

• at* 
MM 

cimmmmmmmmmmmmimmimimtimimmmmi 

RLAltO JAC 0 F. JP.JN. JHA» . 1 . 

I'lNlNElON M < * » 2 A ) t pF£l < 3 * 2 A ) tR'J < 6 » 2 A > »T £ ( 6 * 3 A i 
M MENU ON ueii\il).M(lU!.lt.i| (} 7 t ; » I < 6 » 2 A ' *SM 31*76 J »F < E* 2 A ) 
J»S*.T 3 < 576 > 

COUPON /f Ef'3 / HTtflifd .< i «RU VU,ifl .AU 
/Flfi/ K jhli.'i f>' jr a x, • t jpi N * i jnr.y 

✓ SPACE/ If » NC- 1 f t E C E . ? ) . PSpf < E . 9 > * PEH • f » 9 I • W.M K V J 
✓TITLE / lITCf <30).PAX«.LI1C>» .OELf »«NdlE 

: r_ /?✓ 

:/ 0 / 


NAPE 

1>».I<£ 

t-M'f 

t'U-L 

I'Yl'f 

y 

FAC 

n 

Il- 

ls 

fix 

JRC0L 

hAYG 

NuX 

MYFE 

NTYfE 

l 


TYtf 

r.f At *i 
Hnic 
Uaiis 
MAL*& 
fcl AttP 
W AC » z 
Ifni ci i. 
1 N T E6I.N 
INI I Cl f 
JNIEOI F 

JfUEOlk 

El ALII 

jintCKfi 

I til I Of l- 

JtUtiN 

ItiTCCCR 

REALtO 


SU1R0U1 INIS 


CALIEP: 
LflAT 
CEHf J 


I < y 7 ✓ I- lllf) 

»• t X » / 1; tf£JJ 

r n )/i <eIai 

0 :.T <1 £1 . 

NPl-At K'ln FORCE? 

AXUYKIU.lt IpUOF A1H.N fAll'f 
I'D loot It: Ml 
P 0 l 0 'J! 

I'D l DDE 1 Ml X 
lURREKl ttlriM irri 

».iPiii£Ul' iiirntn im 

I'CllfPlcMl OC JAIDIIAH hAlkl? 

NAYlhOh INC ONE S^Ak C 

RUK'LtUO Ki'Mtf t>r IMERKATltif M'ltnt 
Clif M. til NATIUAI f«»l 
cr.iutiuii hAXttlAC TUI 
MUIU1KJ C MAI If X 


KEH'RNl lHf. t HAIM) 

INITIALIZE THE CNAC E tUNCTll/NS 

MAlR)* HUll UtlCATIDK- 

MATRIX TRAfiSFPSi HUl 7 It 1 1 CAT 1 ON 


hm 
c 


COhfU'N 
COfiPDN 
COUPON 

§*w "A 

Jf11N-9.99li*20 
JMAX--JMIN 

SET OF ELEMENT INf OMIAT lOh Ull 


• OR. *6 .*E. NCX J CALL C6HF 1 UEXtNGXI 


0056 

0059 

$ 82 ? 
00*3 
00*3 
OOtA 
0065 
0064 
067 
066 
069 
0076 

0071 

0072 

0073 
007A 
0075 
O07« 


Emm 

c 

22 


NP-NPGRECt IPX 
Ifllt .ME . JES 
«P 2 >ND*MP 
PC 1 11 -ltNP? 

SMI1 »-0.0 

ZERO OUT POI'T FORCE VECTOR •*** 


t-lilE 

i 


•^?°Sf F rFE, 


CALL CNAT INTTPt ) 


<N7T 

00 10 Il«l»N6 
p«o.o 
Bxi«f -o.o 
PXHE-0.0 

onr-o.o 

DTDE-0,0 

l-MX<l 2 U?EC<I 2 »Il> 

IiXPf -X(12)»I'SHf tI2.IlM0Xl.P 
0XPE-X(12 1*0500(12 I II Hrxur 
pt0f'*t<l'.-)t0ri-Hl2.llMI | TH 
OTOE»TtI2)*PEOE(12»ll I4PY0E 
00 6 12 - 1 i IE 

OrEE<»»I 2 l»P 6 PMJ 2 ilJ)*Pri>E-OSI»E(l 2 ,X*)»PYOF 
““ .< 2 rl 2 )»pXpC*l'SI'E(p»ll)-P*PE»PSI'M 12 iIl > 


5 III 


P-PXPP*OTPI -PXPE*PrP» 


34 


35 
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CC76 
0079 
OOFC 

%m 

0083 
C0E4 

coll 

00P4 

oc-e? 

coce 

©0P9 

OPSO 

00V1 

0091’ 

<077 

0094 

0095 
0 0?4 
0097 
009 E 
©099 
0100 

mi 

0103 
0109 
OKI 

0104 
0107 
OlOfe 
0109 
out 
©111 
0112 
0113 
Oil* 

Cl 13 

one 

0117 

one 

0119 

©120 

0121 

Sis:! 

MSI 
012 1 
0127 
012E 

Sis' 

0131 

0132 

0133 

0134 
0133 
Cl 34 
0137 
0136 

0139 

0140 

mi 

0143 

0144 

%\\l 

0147 

0146 

m 

OlM 

lit! 

0155 

0154 

0157 

0158 
0139 
0100 

S I 41 
142 


£***» 

c 


If I .•.M OP ,|1. 

ir <jacop .c«i 

Ifl.'ACOP .LI • 


JHlN. JNlMJAffU- 
J MAXI JhAX».IAt01 
0.0 1 SI Of 4 


CALCULATE 1 HE *h' MATRIX •»»»> 

VLM/AI1 < 11 1/JAL04 

12*1. It 
61 < 1 » 1 3 ) *0 . 0 
n<?. 131*0.0 

61 (5.131*0.0 

t l < 4 » 1 1 > * 0 . © 

< i » 13 '*!•? ti < i . i :• • 

K’l 131*0.0 
6 ( 3 . 1 3 > * 0 . © 
r <4 » i3i*i>rtt '2.i;* 

4 ( z > . l ?. i * o . o 
lu»!3>*o.o 
r<i .i3}*o.c 
F<2. 131*0.0 

t (5.131*0.0 
F (4.131*0.0 


in (iTiii-o.o 


11 <2.131*0. 
61(3.131*0 


n ( 4 . | 3 > *0 • 0 
n (5.131*0.0 
41 (4.131*0.0 
1(1.131*0.0 
4(2. 13)*lif £1(2.12) 
1(3.131*0.0 
r<4»131*DfEE (1,12) 
4(5.131*0.0 
4(4. 131*0.0 
T<1. 131*0.0 
f <2. 131-0. 0 


£(3.131*0.0 
f <4.131*0.0 
P *5.121*0.0 


|<5« 

|(4> 

Ml. 

<2 


f (4.131*0.0 
I3«1 34 J 

41 (i.i2>»m<i2.in»A' 

61 (2.i2>*FEt(I?.31 1*41 
61 <3.I3)*FEE<12.11 )*#»t 
61(4.131*0.0 

IVMIUM 

6(1.131*0.0 
»(2.13)*0.0 

mm 

31*l>FEt (1.121 
3)*PF£E(2.J2> 

31*0. C 
31*0.0 
F (3.131*0.0 
F (4.131*0.0 

FSS:tti 3 f§U 3 :l!l 

£onJ?NUE 

mh) , rtME „„ 
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IF (JE.HE.O) RETURN 
JE *1 

CONTINUE 
RETURN 
END 
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tmt 


6 

10 
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I 
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©030 

©031 
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©034 

©O3‘o 

0034 

C0J7 

oc!t 

C03* 

0040 

©04) 

0042 

0043 
004 4 
©045 
004© 
©04? 
004C 

0&5C 

©051 

0053 

SE?i 

0055 

0056 

te 

0059 
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A < | * J ) “•© • 0 
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K! I«<1 ./lAHhiF.lF.f 1 )>-«!. /AiKil) 

S IF (rVIWIYIEM 6 0. 40.60 
40 ETA*©.© 

00 IP 6 

60 ALElAXHtl CINIIPEIKCF IHlYf E > »YVIN1tf £>>/!?.» YKIFJT YFE 1 > 
E7A'<l./TANH(AiriA>>-<1 ./ALMA) 

6 r.HA1-<H/2. timSUMJIKim ilCf INTYf E)»XV(N1Y?E > > 4 (£1 AtRG ( Hi YF £ .*» 

*Cr(NTYFE)*YV(KlYlt))l 
10 WI. JTE <6»1©0> KHA1 
101.1 F PFhAT t 4 x . ' F.N61 = ' >2K»G10.5> 

A <5 » 5 ) *KHM 

r l5.5)*>.ViKYYFl >/NV 
(5.6)*YV!Hirri l/NV 
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fc< 6 1 6 >*XV' Hi YF E > /HV 
CALL BHrFt'(A»E»ll»©i6.6) 

CALL Elf Flu*. *1 .MV. 6. 6. 6) 

70 CONTINUE 

Full FLAKE STRAIN »»*» 

C X«EYtN7YIE>/<l .4PMWTYFE1 >/< J .-7 .trMHIVFE) > 

Kl*X6U.-F-MN1rrt} 1 
X7*X*rMH1YF E ) 
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200 rORHATI4X» 'A'.r' »2X»Ci©.5) 

M1>I)*X1 
f. ( 2 » J > » X ? 


.trMHlYFEl? 


Fit .2)«X2 

Mf )»X»U.-2.«FRtNTYPE>>t.5 
K(5.5>»Xl'N1YFE)4UFWI5.1) 
M5.©)*UFU(5.6> 
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RCPV<6»A}»R0<NTYPE>»CFtNWE>»rV<NTYPE> 

RETURN 

CNIi 


Appendix III 
INPUT DATA FORMAT 


Card 

type 

Columns 

Format 

Description 

1 

1-80 

20A4 

Title to be printed at beginning of 
output . 

2 

1-5 

15 

(Number of elements in this problem. 


5-10 

15 

Number of nodes in this problem. 

3 

1-10 

G10.4 

Young's modulus for material type 1. 


11-20 

• 

G10.4 

• 

Yount's modulus for material type 2. 


• 

51-60 

« 

G10.4 

• • • 

• • • 

Young's modulus for material type 6. 

4 

1-10 

G10.4 

Poisson's ratio for material type 1. 


11-20 

• 

G10.4 

• 

Poisson's ratio for material type 2. 

• • • 


• 

51-60 

• 

G10.4 

• • • 

• • • 

Poisson's ratio for material type 6. 

5 

1-10 

G10.4 

Expansion coefficient for material 
type 1. 


11-20 

G10.4 

Expansion coefficient for material 


• 

• 

tvno O 
w J r ^ 

• • • 


• 

51-60 

• 

G10.4 

■ • • 

• • • 

Expansion coefficient for material 

! 



type 6. 

i 6 

1-10 

G10.4 

Conductivity in x-direction for 
material type 1. 


11-20 

G10.4 

Conductivity in x-direction for 


• 

• 

material type 2. 

• • • 


• 

• 

51-60 

• 

• 

G10.4 

• • • 

• • • 

Conductivity in x-direction for 


37 


38 


material type 6. 


7 

1-10 

G10.4 

Conductivity in y-direction for 
material type 1. 


11-20 

G10.4 

Conductivity in y-direction for 


• 

• 

material type 2. 

• • • 


• 

• 

51-60 

• 

G10.4 

• • • 

• • • 

Conductivity in y-direction for 
material type 6. 

8 

1-10 

G10.4 

Density for material type 1. 


11-20 

• 

G10.4 

• 

Density for material type 2. 

• • • 


• 

51-60 

• 

G10.4 

• • • 

• • • 

Density for material type 2. 

9 

1-10 

G10.4 

Specific heat for material type 1. 


11-20 

• 

G10.4 

• 

Specific heat for material type 2. 

• • • 


• 

51-60 

• 

G10.4 

• ♦ • 

• • • 

Specific heat for material type 6. 

10 

1-10 

G10.4 

Velocity in x-direction for material 
type 1. 


11-20 

G10.4 

Velocity in x-direction for material 
type 2. 


• 

• 

• • • 


• 

• 

• • • 


• 

• 

• • • 


51-60 

G10.4 

Velocity in x-direction for material 
type 6. 

1 i 

1-10 

n -i r\ / 

Velocity in y— direction for material 
type 1. 


11-21 

G10.4 

Velocity in y-direction for material 


• 

• 

type 2. 

• • • 


• 

51-60 

• 

G10.4 

• • • 

• • • 

Velocity in y-direction for material 
type 6. 

12 

1-10 

G10.4 

Element height for material type 1. 


11-21 

G10.4 

Element height for material type 2. 


39 


• 

51-60 

• 

G10.4 

• • • 

Element heigh for material type 6. 

1-5 

15 

Node number 

9-9 

11 

Displacement fixity in x-direction at 
this node 

= 0 applied traction 

= 1 applied displacement 

10-10 

11 

Displacement fixity in y-direction at 
this node. 

11-11 

11 

Temperature fixity at this node 
= 0 given flux 
=1 no temperature 

12-21 

G10.4 

x coordinate of node. 

22-31 

G10.4 

y coordinate of node 

32-41 

G10.4 

Force at node in x-direction 

42-51 

G10.4 

Force at node in y-direction 

52-61 

G10.4 

Flux or temperature at node 

1-5 

15 

Element number 

6-10 

15 

Node number for first node on element 

11-15 

• 

15 

• 

Node number for second node on element 

• • • 

• 

41-50 

• 

15 

• • • 

• • • 

Node number for eigth node on element 

57-57 

11 

Material type for element 

64-64 

11 

Element type 


Data cards can be omitted for nodes equispaced between node IT 


and Nj as long as the data for nodes and are included. 
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Example of input data 


r-8E Ffc'iH SEAL 

100 340 i 

.200CL+070.2000E +070 . E00PE+ 07C . 1 775E +0B0 . 1 560E + Of 0 . OOOOE+OO 
.2500 0.2600 0.270C 0.2800 0.3000 0. OOOOE+OC 

.«t3C-L-0: : 0.?7GOr-eOO.e3tOE-OE0.9520£*OM , .i270l-P*0.eOOCi'+C-'. 
.tt iOE-030, J630K-C 7<>. 2C85E-O20. 2552E- C20 . 4 2’2<>E -0 1C. OOPCE+ 

. r 3 2 H -r*3 0 - J 1 3 OF -C .'<< . 7C-yM - (TO . 2552E * 02 3 . 4 2 * ft - C 1 (• . 0 PCC £ i 0 0 
o.it c.pci. o.25< o.?r 


,1:2. 

, j.M OF 
/ . C 

i 


•• .lejo i - . j i j PL -•320 . : ee oe-.-io . 1 55 :>e - c 2 

«.n t> . 0 

000 .ooootfot' 4.850 
C0C'.O<0Cl + ’.'> 4.S8 7 
000 . OOOOF + JO 
00 1.046 


3 ' 

j r 
u 
1 4 

•S * 

* • 

it 

i" 

IF 

1‘ 

20 

22 

32 

:<< 

•»r t 

2# 

•>■> 

2i; 

:■<, 

3: 

23 

i ■ 

37 

3 3 

35 

36 

■3 7 

1 st 

39 

4 0 

4: 

4 2 

4 ;" 

4 A 

45 

44 


00 1 . 05 * 

00 1.86r, 

C r ’ . c “t 

00 : .‘'to 

CO 2 .?ot- 
, ■>? 

66 2.42V 
00 2.4*. j 
00 3.4'-. 
or 2. *71 

00 4.027 
OC 4. ’‘I 
00 4.418 
00 4.440 
00 4.483 
00 4.741 
CO 4.8fO 
0<) 4.880 
C( 4.93C 
OOP. OOOOE+OC 4.925 

ooc.oor OF +00 4 .V 40 

0 • 1.045 


00 1.886 
00 2.002 
0( 2.632 
00 3.483 
CO 3.493 
00 4.047 
00 4,459 
OC 4.472 
CO 4.770 
00 4.®25 
00 4.940 ... 

OOP.OCOOE+OO 4.955 
OCO.OorcE+OO 4.970 
00 1.077 


4.0 e . 0 

G.CCOOt+PQC .CCOCE+OO 
( 1 . oopoe + it 0 .obolti 
o.ooooEit oc.occpe+cc 

0 . OOOCEi POP . OvCCF +00 
<• . OOOOE +000 . C OOOE + 00 

o.oooof+oo;-. opope+cc 
0 . 0 0 00E + 0 to. OOOJf to 0 
4 . (>r-O 0 c + C'O*' . 6 >C 0 CE + C"' 
O.COOOE+OOC . 20006 . 4 c-i 

. OC-JOE + CCC. OCCOE+CO 
0 . 000 PE + 000. C'OOOE + C'O 

o.cocc.ncoo.ocoof+cc 

O.OOOOt+C'OO.CCOOE-fOO 
O.OOPPE+OOO.COOCE+OO 
C . OOOCE + OOP . OOOCE + 00 

O. POOOE+OOP.OOOCE+OO 

P. OOOOE+OOO. OOOOE+OO 
0. OOOOE+OOO. OOOCE+ CO 
c.propE+ocp.ppooi+oo 
O.OOOCt+OPO.OOODE+OJ 

C.tOOOE+POP. OOOOE+OOO. OOOOE+OO 

O . OO0OE + OCC . OOOOE + OOO. OOOOE + 00 

P . OPPOF + OJO . POOCF + POP . OCOCE + OO 
>:■ . cpoOE + poo. 0000 E + 00 
0. OOOOE+OOO. POOOE+DO 
C . OOOOE+OOO .OOOOE+CO 
t*. OOOOE + OOO. OOCOE + OC 
O.POPPEtOOP.PPOOE+OO 
0. OOOOE+OOO. ODDPE+OO 
0. OOOOE+OOO. OOOOE+OO 
O.OOOOf+OOO.CCOOE+OC 


4.910 
4.6E2 
4.741 
4 . 3 - P 3 
4.43P 
4.445 
? . 9 ? 2 
4 . 02 :.' 
:-.4?9 

2.453 

3.477 

2.760 

•» — (5 * 

1*966 
1.978 
J .990 
1 .046 

1 . 05 V 


4.770 

4.459 

4.472 

4.047 

3.483 

3.492 

2.61? 

1.996 

2.002 

1.063 


g.pOOOE+OOO.OOOOE+OO 


00 2.008 
00 2.014 
OC 2.829 
00 3.504 
00 3.514 
OC 4.071 


47 

CO 4.466 

2. COE 

48 

00 4.499 

2.014 

4 V 

00 4.789 

5.072 

»/l 

00 4.855 C 

.OOoOt 

f. - 

00 4 . 8 7 C C 

.00 JOE 

K •• 

OOP . OOOOE+OO 

4.985 


000. POME + 00 

5.00C 

5 « 

00 3 . 

4 .828 

55 

00 2.021 

4.513 

56 

00 2-027 

4.527 

• ,% 

00 7.64/, 

4.096 

58 

OC 3.525 

% . r .v* 


00 2.536 

3.536 


_ OOOOE+OOO. OOOOE+CO 
0. OOOOE+OOO. OOCOE+OO 
0. OOOOE+OOO. OOOOE+OO 
O.COOOE+OOC. OOOOE+OOO. OOOOE+OO 
0. OOOOE + OOO. OOC’CE + OOO. OCOOE + OO 
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